An Analytic Theory for the Orbits of Circumbinary Planets 

Gene C. K. Leung 1 and Man Hoi Lee 1,2 

ABSTRACT 

Three transiting circumbinary planets (Kepler-16 b, Kepler-34 b, and Kepler-35 b) 
have recently been discovered from photometric data taken by the Kepler spacecraft. 
Their orbits are significantly non-Keplerian because of the large secondary-to-primary 
mass ratio and orbital eccentricity of the binaries, as well as the proximity of the planets 
to the binaries. We present an analytic theory, with the planet treated as a test particle, 
which shows that the planetary motion can be represented by the superposition of the 
circular motion of a guiding center, the forced oscillations due to the non-axisymmetric 
components of the binary's potential, the epicyclic motion, and the vertical motion. In 
this analytic theory, the periapse and ascending node of the planet precess at nearly 
equal rates in opposite directions. The largest forced oscillation term corresponds to a 
forced eccentricity (which is an explicit function of the parameters of the binary and 
of the guiding center radius of the planet), and the amplitude of the epicyclic motion 
(which is a free parameter of the theory) is the free eccentricity. Comparisons with 
direct numerical orbit integrations show that this analytic theory gives an accurate 
description of the planetary motion for all three Kepler systems. We find that all three 
Kepler circumbinary planets have nonzero free eccentricities. 



INTRODUCTION 



Doyle et al.1 (|201ll ) have recently discovered the first transiting circumbinary planet, Kepler- 
16 b, from photometric data tak en by the Kepler sp acecraft. The Saturn-mass planet orbits a pair 
of stars of O.69M and 0.20M Q . I Welsh et al.1 (12012] ) subsequently announced the discovery of two 
more circumbinary planets: Kepler-34 b and Kepler-35 b. Kepler-34 b is 0.22Mj (where Mj is the 
mass of Jupiter) and orbits two Sun-like stars, while Kepler-35 b is 0.13Mj and orbits a pair of 
smaller stars (O.89M0 and 0.81Mq). For all three systems, the orbits of the binary a nd planet are 



aligne d to within 2°. From the observed rate of circumbinary planets in their sample, IWelsh et al 
(|2012j ) estimated that more than ~ 1% of close binary stars have giant planets on nearly coplanar 
orbits. 

Variations in eclipse times and transit durations, combined with radial velocity measurements, 
allow precise measurements of both physical and orbital parameters for all three systems. Tabled] 
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shows the best-fit osculating Keplerian orbital parameters provided by J. A. Carter (2012, private 
communication). They differ slightly from those published in Table 1 of I Welsh et al.l (|2012l ). as the 
values in that table are the medians of the cumulative distribution of the marginalized posteriors 
for each parameter, and they are osculating parameters at a different epoch. 

Plots of the evoluti on of the osc u lating Keplerian orbital elements of the pla nets in the Support - 
ing Online Material of iDovle et al.l (|201ll ) and Supplementary Information of IWelsh et al.l (|2012l ) 
show significant variations on both orbital and secular timescales, with the eccentricity changing 
from nearly zero to 0.1 in the case of Kepler-16 b, and the precession period of the orbit is as short 
as ~ 60 orbital periods in the case of Kepler-35 b (see figures below for more details). The non- 
trivial departures from unperturbed Keplerian orbits for these circumbinary planets are due to the 
large secondary-to-primary stellar mass ratio (mg/m^ = 0.29-0.97), the large orbital eccentricity 
of the binary eAB = 0.14-0.52, and the proximity of the planet to the binary (orbital period ratio 
Pb/PAB = 5.6-10.4). 

In £}2] we present an analytic theory for the orbit of a circumbinary planet that is valid in the 
limit that the planet has negligible mass and can be treated as a test particle. iLee Peald ((2006) 
have previously developed an analytic theory for the orbits of the small satellites of the Pluto- 
Charon system, assuming that the orbit of Charon relative to Pluto is circular. We generalize that 
theory to the case of an eccentric binary orbit. In $3] we present the results of direct numerical orbit 
integrations and compare them with the analytic theory. In §4] we discuss the limitations of the 
analytic theory and a simple modification that can significantly improve the analytic predictions 
when the orbital eccentricity of the binary is large. Our conclusions are summarized in Sj5j 



ANALYTIC THEORY 



In this section we follow a similar approach as lLee Peald (|2006l ) to develop an analytic theory 



for the orbit of a circumbinary planet. We extend their theory to first order in the eccentricity of 
the orbit of the binary. We assume that the planet has negligible mass and can be treated as a test 
particle. Then the orbit of the secondary (hereafter B) relative to the primary (hereafter A) is an 
elliptical Keplerian orbit with eccentricity cab an d semimajor axis ciab , and the distance between 
A and B is tab = oab(1 — e AB )/(l + eAB cos f B ), where f B is the true anomaly of B. We adopt 
a cylindrical coordinate system with the origin at the center of mass of the binary and the x-y 
plane being the orbital plane of the binary. The positions of B and A are r B = (Rb, 4>b, 0) and 
r A = {Ra,4>b + vr,0), respectively, where R B = f"AB m A/ ( m A + m B ), Ra = r AB m B/ i m A + mn), 
m,A is the mass of A, and m B is the mass of B. 
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2.1. Potential 



The gravitational potential at r = (i?, 0, z) due to the binary is 



$(r) 



(1) 



Since the orbit of the planet is nearly coplanar with that of the binary, we expand l/\r — r B \ in 
powers of z: 

1 "V", (2) 



where 



\r-r B \ ( p 2 + z 2)i/2 p 2p 3 
p = [R 2 + R 2 B - 2RR B cos (<j) - <j) B ) 



1 1/2 



(3) 



The in verse powers of p can be expressed as cosine series using equation (6.62) of lMurrav Dermott 
(|l99gh to give 
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where S^o is the Kronecker delta and b^(R B /R) are the Laplace coefficients. Similarly. 

1 ^E(-i)*(2-fc,) r - ~ "~ 1/z ' 2 
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Then the potential can be written as 

oo 

*(r) = 
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where 
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To first order in e AB , 



= !b + ~ Mb + sin M B + ro B , 



(7) 



(8) 



where w B and M# = n AB t + yj^s are, respectively, the longitude of periapse and mean anomaly 
of B relative A, n AB = [G(m A + ms) / a AB \ l l 2 is the mean motion of the binary, and (f AB is a 
constant. Then 



cos k(4> — 4> B ) ?» cos — Mb — w b) + e^fcfcos 
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To first order in cab, 



Rb/clb ~ 1 - e A B cos M B , 



(10) 



and 



b k {j+1)/2 (R B /R) « ^ 7+ i)/ 2 ( q; b) - eABa B D&( - +1 y 2 (a B )cosM B , 



J (j+l)/2\ sh B/^) ~ »(j+l)/2V«B; _ c^B« J B-Ly^(j + i)/ 2 V«^;^ UBi,/J -B; (H) 

where a# = aABmi/(m/i + tub), «b = as/R, and D = d/da. Equation (fTTj) gives rise to terms 
involving 



&AB cos Ms cos k{4> — 4>b) 



eAB 



[cos(fc(0 - <f> B ) + M B ) + cos(fe(0 - <f> B ) - Mb)] , (12) 



which can be expressed as 

cab cos(k(<j) - 4>b) ± M B ) ~ eAB cos(k(4> - w B ) - {k T ^)M B ) 
Terms involving Ra/R can be found in a similar manner. 
After re-grouping the terms, 
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and 
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The terms in equation (1141) multiplied by eA R are new compared to the potential due to a binary on 
circular orbit derived bv lLee Peald (|2006l ). As in the circular binary orbit case, the axisymmetric 
$j00 components of the potential are identical to those due to two rings: one of mass tua and 
radius a a and another of mass m B and radius a B . 



2.2. Equations of Motion and Solution 

The equations of motion in cylindrical coordinates are 

12 d$ 



R-Rtf 
R4> + 2R4> 



dR' 
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R~d0' 
d$ 

dz 



(17) 



As in lLee &z Peald (j200Q ). we approximate the orbit of the circumbinary planet as small deviations 
from the circular motion of a guiding center in the x-y plane: 

R = Ro + Ri{t), 

<t> = Mt) + Mt), (is) 

z = z 1 (t), 

where the constant Ro is the radius of the guiding center, \Ri/Rq\ <C 1, \<pi \ <C 1, and |zi/i?o| ^ 1- 

Substituting equations (TLT1) and (fT8l) into equation (fT7|) . the only nontrivial equation at the 
zeroth order is 

~d®ooo~ 



. 2 

r 4>o 



dR 



Ro 



which describes the circular motion of the guiding center. Its solution is 

4>o (t) = n t + ifo, 
where <po is a constant and the mean motion no is given by 
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In the above equation uk = [G(m A + m B)/^o]^ 2 is the Keplerian mean motion at Rq, and a a 
and a B are evaluated at R = Rq. 
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At the first order, the equations of motion are 
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where n = (R~ 1 d<&Qoo/ dR) 1 / 2 is the mean motion at R, and the quantities in the square brackets 
with the subscript Ro are evaluated at R = Rq. Equation (j23|) can be integrated to give <pi, which 
can then be substituted into equation ([22]) to yield 
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where the epicyclic frequency kq is given by 
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Equation (j25j) is the equation of motion of a simple harmonic oscillator of natural frequency kq 
that is driven at frequencies n AB , k\no — n AB \, and \kno — (k ± l)nAB|, and its solution gives 

f oo 
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where ef ree and are arbitrary constants and 
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We can then solve equation (|23|) to give 
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Dr. 
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The motion in R and <f> is the superposition of the circular motion of the guiding center at Rq at 
frequency no, the epicyclic motion represented by the free eccentricity ef ree at frequency kq, and the 
forced oscillations of fractional radial amplitudes Co, C° and at frequencies uab, k\riQ — uab\ 
and \kno — (k ± ^n^l, respectively. Note that Co, C°, or C^ 1 is singular if Kq — n 2 AB , kn — Iuab, 
or Kq — (knQ — Iuab) 2 = 0, where I = k or k rb 1. The second and third com binations of frequencies 
corre spond to the corotation and Lindblad resonances, respectively (e.g., iGoldreich Tremaine 
19801 ). None of these resonances could be encountered if the planet is further away than the 3:1 
mean-motion resonance with the binary. 



(35) 



For the motion in z, the solution to equation (|24j) is 

Z = Z\ = Rokrce COs(u t + (), 

where if ree and £ are arbitrary constants and the vertical frequency z^o is given by 
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Thus the motion in z decouples from that in R and 
the free inclination if ree at the vertical frequency uq. 



{m A + m B ) 
and has only free oscillations represented by 



As we shall see, circumbinary planets have vq > hq > > kq. So the azimuthal period 2-7r/no 
is shorter than the Keplerian orbital period 2n/riK, the periapse precesses in the prograde direction 
with the period 2-7r/|no — kq\, and the ascending node precesses in the retrograde direction with the 
period 2tt /\uq — 



As in the circular binary orbit theory of lLee &: Peald (J2006J), the motion is represented by the 
circular motion of the guiding center, the epicyclic motion, the forced oscillations and the vertical 
motion. The expressions for the mean motion no (eq. [H]), the epicyclic frequency hq (eq. |26j). 
and the vertical frequency vq (eq. [36]) involve the axisymmetric $000 an d ^200 components of the 
potential and are identical to those in the circular binary orbit case (there are, however, corrections 
at the second order in eAB] see The motion in z is identical to the circular binary orbit case. 
The forced oscillations are composed of both terms identical to those in the circular binary orbit 
theory (C° at frequencies k\riQ — uab I) and new terms (Co at frequency uab and at frequencies 
\kno — (k ± 1)iiab\)- Note that the new terms Co and C k are proportional to eAB and that one 
of these terms, Cj~ , has frequency uq and can be identified as the forced eccentricity. The forced 
longitude of periapse is aligned with the binary's periapse because Cf > 0. 
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If we expand the analytic expressions in powers of clab/Ro (note that clab/Ro % 0.32 for the 
three Kepler systems) and retain only the lowest order term, we find that the precession rates 

2 



n - k _ n - up ^ 3 m A m B ( oab\ 
n K n K 4 (m A + m B ) 2 \ Rq J 

the forced eccentricity 
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We can see from equations ([29]) and ([30]) that C° and Cj*r involve ^ofcO) &oki, and their derivatives 
with respect to R. According to equations f|15|) and f 1 1 6 [) . these terms would be exactly zero if k is 
odd and txia = m B - Equations (|38|) . ([41 p . and (j42p show that the odd terms are proportional to 
{mA — ni B )/(m,A + m B ) at the lowest order in clab/Rq and could be small if uia ~ tti b . 



COMPARISONS WITH NUMERICAL ORBIT INTEGRATIONS 



For the numerical orbit integrations, we use Jacobi coordinates (where the position of the 
secondary B is relative to the primary A and the position of the planet is relative to the center of 
mass of the binary), with the invariable plane as the reference plane. This coordinate system reduces 
to that used in £j2] when the mass of the planet is negligibl e. We perform dire c t num erical orbit 
integrations of the Kepler-16, 34, and 35 syste ms, using the lWisdom fc Holmanl (jl99ll ) symplectic 
integrator with the modification described in iLee &: Peald ([20031 ). The modification allows the 
integration of circumbinary planets without an excessively small timestep, and we use a timestep 
of 0.1 day. We generate the initial positions and velocities of the binary and planet by using the 
best-fit osculating orbital parameters in Table [TJ 

For comparison with the analytic theory, we need to evaluate nx , no, ko, vq, and the fractional 
amplitudes Co, C°, and C^ at a guiding center radius Ro- For each system, we adopt the average of 
the maximum and minimum values of the cylindrical radius R of the planet's orbit in the numerical 
orbit integration over many precession cycles as Rq- The adopted Ro and the numerical values of 
n Ki n o/nK, kq/tlk, vo/tik, Co, Cj?, and (k = 1, 2, and 3) are listed in Table [2] 
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3.1. Kepler-16 



In Figure Q] we plot the evolution of the osculating Keplerian orbital elements of the planet 
Kepler-16 b over 100 years from the numerical orbit integration. The eccentricity eb shows variations 
on both orbital and apsidal precession timescales. The longitude of periapse Wb changes rapidly 
when eft is nearly zero, but the long-term trend is prograde precession with a period of 48.6 years. 
The precession of the longitude of ascending node Qb is retrograde with a period of 41.0 years, 
and the inclination % shows small oscillations around a constant value, with two oscillations per 
nodal precession period. Using nx, no/nx, k-q/uk, and vq/uk from Tabled the analytic theory 
predicts that the apsidal and nodal precessions are at nearly equal rates in opposite directions, with 
the prograde apsidal precession having a period of 42.2 years and the retrograde nodal precession 
having a period of 42.8 years. These are in good agreement with the numerical results but slightly 
shorter for the apsidal precession and longer for the nodal precession. 

In the bottom panels of Figure [2] we plot the variations in the orbital radius Rb of the planet 
Kepler-16 b over 5.4 years and 100 years. Significant periodic variations in the amplitude of the 
oscillations in Rb are observed from the 100-year plot. The variations have a period of 48.6 years, 
which is equal to the period of apsidal precession. The variations are the result of the superposition 
of the epicyclic motion at frequency kq with amplitude ef ree and the forced oscillation at frequency 
no with amplitude Cj~. Without any loss of generality, we can assume that both ef ree and Cj~ > 0. 
After some algebraic manipulation using sum and product formulae of trigonometric functions, 
these two terms in equation (|27p can be written as 



where Wf vee = (no — Ko)t + (fo — ift is the free longitude of periapse. A maximum amplitude is 
reached when Wf ree — wb = 2£ir (where i is an integer), in which case the right-hand side of 
equation (|4"3]) becomes ±(C{" + ef ree ) cos{[(no + Ko)t + ipo — wb + ^]/2}. Similarly, a minimum 
amplitude is reached when Wf ree — wb = (21 + l)vr, in which case the right-hand side of equation 
([4*3]) becomes ±(C{~ — ef ree ) sin{[(no + Ko)i + <A) — WB + ift]/2}- Therefore, a maximum (or minimum) 
amplitude is reached every 2n/\no — kq\, which is equal to the apsidal precession period. The small 
minimum amplitude and large variations in the amplitude observed in the 100-year plot indicate 
that ef re e ~ C±- 

The 5.4-year plot in the lower left panel of Figure [2] clearly shows an increase in the amplitude 
of the oscillations in Rb at the initial epoch due to the changing relative phase between the free 
and forced eccentricity terms, as well as higher-frequency forced oscillations. In order to study in 
more detail the forced oscillations and epicyclic motion, we plot in the upper panels of Figure [2] a 



efree COs(/-v0* + 1p) + C 1 COs((po — Wb) 

= efree cos(ko^ + ip) + cos(not + ifo — Wb) 




(n + Kp)t + (po-WB + l/j ~ 

2 

(n + Kp)t + (fo - w B + tp 
2 



(43) 
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transformed orbital radius denned by 



R' = R + Rq^Cq cos M B + J2[ C k cos k{cj)o - M B -zu B ) 



+C+ cos^^o - vj b ) ~(k+ l)M B ) + C~ cos(k{(/> - w B ) - (k - l)M B )\ j , (44) 

with i?o, Co, C£, and from Table [2] and 4>q, Mb, and wb from the numerical integration 
itself (which eliminates phase errors due to small frequency errors and the very slow precession 
of the binary's periapse). It is clear from a comparison between the upper and lower panels of 
Figure [2] that the forced oscillations (including the forced eccentricity term) are sufficiently close 
to those predicted by the analytic theory that they are effectively eliminated in R' b . The largest 
forced oscillation term, other than the forced eccentricity term, is with a fractional amplitude 
of 0.0024 and a period of 27r/(2n — ti^b) = 64.5 days, and the other forced oscillation terms are 
at least a factor of 4 smaller in amplitudes. The free eccentricity can be easily obtained from R' b , 
which shows only periodic epicyclic variation at frequency kq. The forced and free eccentricities of 
Kepler-16 b are = 0.036 and ef ree = 0.030, respectively. 

Having obtained ef ree , which is a free parameter of the analytic theory, we can now directly 
plot the evolution of R according to equation (p7|) of the analytic theory (Fig. [3j), as well as the 
evolution of the osculating Keplerian elements using R from equation ([27]), <p from equation (I3ip . 
and their time derivatives (Fig. |4|). They are in excellent agreement with the numerical orbit 
integration shown in Figures Q] and EJ except for the faster periapse precession (and hence faster 
periodic variations in the radial oscillation amplitude), the slower nodal precession, and the lack 
of periodic variations in the inclination. Our analytic theory only gives the free oscillations in the 
vertical direction and cannot explain the periodic variations in the inclination at twice the nodal 
precession rate. For the eccentricity variations, if we ignore the higher-frequency forced oscillations, 
equation (jl3|) shows that the osculating eccentricity reaches a maximum of C{" + ef ree = 0.066 when 
the free longitude of periapse zo{ me is aligned with the forced longitude of periapse (which is equal to 
the longitude of periapse of the binary wb\ see £ j2.2p . and reaches a minimum of \C± — ef ree | = 0.006 
when the free longitude of periapse rof ree is anti-aligned with the forced longitude of periapse. This 
behavior agrees with the usual definitions of the forced a nd free eccentricities and longitudes of 



periapse (see, e.g., Section 7.4 of Murray k, Dermottlll999l ). 



3.2. Kepler-34 

Figures [5] and [6] show the evolution of the osculating Keplerian orbital elements, the orbital 
radius i?& and the transformed orbital radius R' b of the planet Kepler-34 b. The periods of prograde 
apsidal precession and retrograde nodal precession are 62.9 and 67.9 years, respectively, from the 
numerical orbit integration. The analytic theory predicts 91.1 and 91.9 years, respectively, which 
are longer than the numerical results by more than 35%. The main reason for the large discrepancy 
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is that the analytic theory is only accurate to first order in the binary eccentricity cab and Kepler- 
34 has the largest 6ab(= 0.52) among the three systems. We shall derive in 2] simple corrections 
at the second order in cab that significantly improve the analytic predictions of the precession 
periods. 

One might think that the large osculating Keplerian eccentricity (e& ~ 0.2) is due to forcing by 
the eccentric binary. But the nearly identical plots of Rb and R' b show that the variations in are 
primarily due to epicyclic motion and that the forced eccentricity and other forced oscillation 
terms are small. Indeed, Cj~ = 0.0019, which is smaller than that for Kepler-16 b by more than 
an order of magnitude, and the next largest forced oscillation term is = 0.00068 (see Table [2]). 
The forced eccentricity C{~ , as well as all C£ and terms with odd k, are small because of the 
nearly equal masses of the binary components of Kepler-34 (m^/m^ = 0.97; see discussion in the 
last paragraph of §2.2|) . We find from the variations in R' b that the free eccentricity ef ree = 0.204. 

3.3. Kepler-35 

For Kepler-35 which has binary eccentricity {eAB = 0.14) similar to Kepler-16, the numeri- 
cal and analytic apsidal and nodal precession periods of the planet are in good agreement. The 
numerically determined apsidal and nodal precession periods are 21.7 and 20.2 years, respectively 
(see Fig. [7J, and the analytic ones are 20.4 and 20.8 years, respectively. 

As in the case of Kepler-34, the binary components have nearly equal masses {tub I'm A = 
0.91) and the forced eccentricity of the planet C± = 0.0025 is small. However, because the free 
eccentricity is much smaller than that for Kepler-34 b and comparable to that for Kepler-16 b, 
moderate variations in the amplitude of oscillations in Rf, with the same period as the apsidal 
precession are clearly observed in the 100-year plot in the lower right panel of Figure [HJ The 
4-year plot in the lower left panel of Figure [8] also shows the effects of higher-frequency forced 
oscillations terms. As for Kepler-16 b and 34 b, is the largest forced oscillation term after 
the forced eccentricity term C7 (see Table [2]) . The forced oscillations are sufficiently close to those 
predicted by the analytic theory that the transformed orbital radius R' b shows only periodic epicyclic 
variation at frequency kq (see upper panels of Fig. [8]). The free eccentricity from the variations in 
R' b is e free = 0.038. 

4. DISCUSSION 

The analytic theory developed in $2] is accurate to first order in the binary eccentricity eAB 
and to first order in the deviations R%, <j)i, and z\ of the planetary motion from the circular motion 
of the guiding center. It also treats the planet as a test particle and ignores the gravitational 
effects of the planet on the motion of the binary. From the comparisons with direct numerical orbit 
integrations of the Kepler-16, 34, and 35 systems in $3l we have shown that the analytic theory 



-13- 



gives an accurate description of the planetary motion in all three cases, except for the apsidal and 
nodal precession periods of Kepler-34 b with cab = 0.52. 

It was pointed out in $2] that the expressions for the mean motion no (eq. |21j). the epicyclic 
frequency kq (eq. [26]), and the vertical frequency uq (eq. [36]) involve the axisymmetric $000 
and $200 components of the potential and that the axisymmetric components of the potential are 
identical to those due to two rings: one of mass tua and radius a a an d another of mass m» an d 



radius as- If we expand Rb/cib to higher orders in cab (see eq. [2.81] of lMurrav k, Dermottlll999l ) 



D g2 3 e 3 

_ ^ = l _ eAB cos M B + -^(1 - cos 2M B ) + — ^(cos M B - cos 3M B ) + • • • , (45) 
<iB 2 8 

which means that the time-averaged Rb/clb = 1 + e \ B /2 and that it might be more appropriate 
for the axisymmetric components of the potential to be due to two rings of radii a^(l + e 2 AB /2) 
and a#(l + e AB /2). This is achieved if we selectively include just the e AB /2 term beyond first 
order in cab and use Rb/clb ~ 1 — cab cos Mb + c\ B /2 (and similarly for Ra/cia)- Then the 
only modifications to the analytic theory in SJ2]are that ^(j +1 )/ 2 ( a ^) ^ s re pl ace d by b^ +1 y 2 [aA{l + 
e As/ 2 )]' and ^\j+i)/2^ aB ) b y tf j+1 y 2 [aB(l + e As/ 2 )]' in equation (Jl5j) for ® jk0 , and similarly 
for Db k ^ +l y 2 (a>A) and DbQ +1 y 2 (ctB) in equation (fl~6|) for Qjki- With this simple modification, 
the analytic predictions for the apsidal and nodal precession rates are faster than the unmodified 
values by only 2-3% for Kepler-16 b and 35 b, but by ~ 29% for Kepler-34 b. The increase by 
approximately (1 + e AB /2) 2 can be understood from the (ciab/Ro) 2 scaling of the lowest order 
expression for the precession rates in equation (|37p . The modified analytic precession periods for 
Kepler-34 b are 71.4 years for periapse and 72.1 years for ascending node, which are much closer to 
the numerical results (62.9 and 67.9 years, respectively) than the unmodified values (~ 91 years). 

The 1 + e AB /2 modification also affects the amplitudes of the forced oscillation terms. The 
change in the largest of these, the forced eccentricity C{~ , is small: < 6% even for Kepler-34 b. The 
change in the second largest forced oscillation term, C 2 , is approximately 1 + 5e^ B /6, which is 
only 2% for Kepler-16 b and 35 b but ~ 24% for Kepler-34 b. However, as we saw in §3.21 both C± 
and C 2 are small compared to the free eccentricity and have no noticeable effect on the evolution 
of Rf, for Kepler-34 b. 

The 1 + e AB /2 modification that we have just described is not rigorously correct. We have at- 
tempted to construct an analytic theory that is accurate to 0(e AB ). Preliminary analysis indicates 
that there are no additional corrections at 0(e AB ) for C 2 but that there are additional corrections 
to the precession periods and C-f . The full 0(e AB ) corrections to the precession periods and C± of 
Kepler-16 b and Kepler-35 b, as well as of Kepler-34 b, remain small (less than a few %). For 
Kepler-35 b, the apsidal (nodal) precession period may decrease (increase) by a few years beyond 
the 1 + e AB /2 modification discussed above. The full 0(e AB ) theory also introduces new forced 
oscillation terms with frequencies \kriQ — (k ± 2)tiab\, but they are likely small in amplitudes, as 
they are not observed in the direct numerical orbit integrations. 



The most obvious effects of the gravitational force of the planet on the binary are the precession 
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of the binary's periapse and ascending node. Due to our choice of the invariable plane as the 
reference plane, the longitude of ascending node of the binary must be 180° from, and precesses at 
the same rate as, the longitude of ascending node of the planet. From the direct numerical orbit 
integrations, we find that the apsidal precession rates of the binary are 0.026, 0.0033 and 0.0086 
degrees per year for Kepler-16, 34, and 35, respectively, which are much smaller than those of the 
planet. For the comparisons in ^3j this very slow precession of the binary's periapse is taken into 
account by using wb from the numerical integrations when we plot R' defined in equation (j44]). 

5. CONCLUSIONS 

We have developed an analytic theory to model the motion of the recently discovered cir- 
cumbinary planets: Kepler-16 b, Kepler-34 b, and Kepler-35 b. Their orbits are significantly 
non-Keplerian due to the large secondary-to-primary mass ratio and orbital eccentricity of the bi- 
naries, as well as the proximity of the planets to the binaries. The analytic theory in ^2] shows that 
the motion in R and <f> can be represented by the superposition of the circular motion of a guiding 
center at mean motion no, the epicyclic motion, and the forced oscillations, and that the motion 
in z decouples from that in R and (f> and has only free oscillations. One of the forced oscillation 
terms has frequency no and can be identified as the forced eccentricity, while the epicyclic motion 
at frequency kq can be identified as the free eccentricity. 

Comparisons with direct numerical orbit integrations in £}3] show that the analytic theory (with 
the simple modification in Q gives an accurate description of the planetary motion for all three 
Kepler systems, including the precession of the periapse and ascending node. The analytic theory 
explains the periodic variations in the amplitude of oscillations of the orbital radius (which is most 
obvious for Kepler-16 b and negligible for Kepler-34 b) by the superposition of the epicyclic motion 
and the forced eccentricity oscillation. The amplitude (and osculating eccentricity) variations have 
a period equal to that of the apsidal precession as predicted by the theory. 

The amplitude, C-f , of the forced eccentricity term is an explicit function of the parameters 
of the binary and of the guiding center radius of the planet in the analytic theory. For Kepler- 
16 b, 34 b, and 35 b, C-f = 0.036, 0.0019, and 0.0025, respectively. The free eccentricity, ef ree , 
of the epicyclic motion is a free parameter in the analytic theory and can be obtained from, e.g., 
the variations in the orbital radius of the planet. For Kepler-16 b, 34 b, and 35 b, ef ree = 0.030, 
0.204, and 0.038, respectively. Note that the Kepler-34 system with the largest binary eccentricity 
{cab = 0.52) has the largest ef ree while the other two Kepler systems with comparable cab have 
comparable ef ree . Since the free eccentricity is a free parameter that was set by dynamical processes 
during the formation and/or subsequent evolution of the circumbinary planet, the free eccentricity 
of the three Kepler circumbinary planets (and additional circumbinary planets in the future) should 
provide important clues to these processes. 

While this paper was under review, three more circumbinary planetary systems were an- 
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nounc ed: Kepler-38 and PHI with one planet each and Kepler-47 with two planets (jOrosz et al 



2012al lbl: [ Schwamb et al.ll2012l ). Direct numerical integrations and our analytic theory show that 



(i) Kepler-38 b is similar to Kepler-16 b in having ef ree ~ ~ 0.024 and large variations in 
the amplitude of oscillations in R; (ii) Kepler-47 b is similar to Kepler-35 b in having ef ree larger 
than C{~ ~ 0.004 and moderate variations in the amplitude of oscillations in R; (hi) Kepler-47 c is 
similar to Kepler-34 in having ef ree much larger than C-j" ~ 0.001 and negligible variations in the 
amplitude of oscillations in R; and (iv) PHI has ef ree ~ 0.1 and ~ 0.04. 

It is a pleasure to thank Josh Carter for furnishing the best-fit orbital parameters of the Kepler 
circumbinary planetary systems, Daniel Fabrycky for informative discussions, and the referee for 
helpful comments on the manuscript. C. K. L. also thanks K. H. Chan, X. Tan, and W. Y. Li 
for enlightening discussion. This work was supported in part by Hong Kong RGC grant HKU 
7034/09P. 
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Table 1. Orbital Parameters of Circumbinary Planetary Systems 



Parameter 




Kepler- 16 


Kepler-34 


Kepler-35 


Epoch (BJD) 




2,455,000.0 


2,454,900.0 


2,454,900.0 


GM A (10- 4 AU 3 or 2 ) 




2.0328 


3.1045 


2.6187 


GM B (10~ 4 AU 3 d- 2 ) 




0.5987 


3.0232 


2.3903 


GMb (10 Air d z ) 




9.3119 


6.5822 


3.6839 




Binary star orbit 






Semimajor axis (AU) 




0.22405 


0.22847 


0.17603 


Eccentricity 




0.16048 


52068 


0.14224 


Inclination (deg) 




0.0011 


0.0020 


0.0006 


Argument of periapse 


(deg) 


257.79 


323.86 


338.95 


Long, ascending node 


(deg) 


5.70 


107.45 


107.56 


Mean anomaly (deg) 




129.84 


52.66 


299.31 






Planet orbit 






Semimajor axis (AU) 




0.72042 


1.08617 


0.60497 


Eccentricity 




0.02373 


0.20861 


0.04845 


Inclination (deg) 




0.3083 


1.8590 


1.0714 


Argument of periapse 


(deg) 


304.05 


69.41 


91.17 


Long, ascending node 


(deg) 


185.70 


287.45 


287.56 


Mean anomaly (deg) 




358.85 


17.75 


292.17 



Note. — The orbital parameters are the best-fit osculating Jacobian 
parameters relative to the invariable plane at the listed epoch. 
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Table 2. Parameters of Analytic Theory 



Parameter 


Kepler- 16 


Kepler-34 


Kepler-35 


Ro (AU) 


0.7016 


1.0804 


0.5933 


riK (yr _1 ) 


10.0823 


8.0512 


17.8875 


n /n K 


1.00702 


1.00423 


1.00838 


K /n K 


0.99224 


0.99567 


0.99119 




1.02158 


1.01272 


1.02527 


Co 


0.000159 


0.000085 


0.000131 


C? 


-0.000282 


-6 x 10~ 7 


-0.000020 


c° 2 


-0.000589 


-0.000079 


-0.000533 


cl 


-0.000049 


-1 x 10~ 7 


-0.000003 


ct 


0.000005 


4 x 10~ 8 


3 x 10~ 7 


ct 


-0.000033 


-0.000016 


-0.000028 


ct 


-0.000006 


-4 x 10~ 8 


-4 x 10~ 7 


Ci 


0.035772 


0.001861 


0.002493 


C2 


0.002438 


0.000683 


0.001731 


Cs 


0.000110 


7 x 10~ 7 


0.000007 



20 40 60 80 100 
t (yr) 



20 40 60 80 100 
t (yr) 




Fig. 1. — Evolution of the osculating Keplerian orbital elements (eccentricity e^, inclination %, 
longitude of periapse m^, and longitude of the ascending node Qj,) of the planet Kepler-16 b over 
100 yr from direct numerical orbit integration. The elements are relative to the center of mass of 
the binary, and the reference plane is the invariable plane. 
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Fig. 2. — Variations in the orbital radius Rb (bottom panels) and the transformed orbital radius 
R' b (top panels; eq. [44j ) of the planet Kepler-16 b over several years (left panels) and over 100 yr 
(right panels) from direct numerical orbit integration. 
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Fig. 3. — Variations in the orbital radius Rt, of the planet Kepler-16 b over several years (left panel) 
and over 100 yr (right panel) according to equation (|27|) of the analytic theory with ef ree = 0.030. 




Fig. 4.— 



Same as Fig. [H but according to the analytic theory with ef ree = 0.030. 
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Fig. 6. — Same as Fig. [2j but for the planet Kepler-34 b. 
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Fig. 7. — Same as Fig. [Q but for the planet Kepler-35 b. 
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Fig. 8. — Same as Fig. [2j but for the planet Kepler-35 b. 



